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Abstract. A method is proposed to design the time dependence of the trap frequency 
and achieve in a short time an adiabatic-like (frictionless) evolution of Bose-Einstein 
condensates governed by the Gross-Pitaevskii equation. Different cases depending on 
the effective dimension of the trap and the interaction regimes are considered. 2D traps 
are particularly suitable as the method can be applied without the need to impose any 
additional time-dependent change in the strength of the interatomic interaction or a 
Thomas-Fermi regime as it occurs for ID and 3D traps. 

PACS numbers: 67.85.De, 42.50.-p, 37.10.Vz 
1. Introduction 

In order to manipulate Bose-Einstein condensates for different applications it is 
important to study and control their response to time-dependent changes of the confining 
fields. A natural approach to avoid undesired excitations is to modify the trap 
adiabatically, i.e., very slowly, so that, if the initial state is in the ground state the 
final state will be the ground state as well. However, this may require very long times 
and become impractical. Faster changes are thus a desirable objective but they will 
in general induce excitations and oscillations (inner frictional heating [1]), so that the 
proportion of the ground state in the final state may be small [21 El Hj. These difficulties 
raise the question addressed in this paper: Is it possible to change the trap in a very 
short time, taking the condensate, up to a global phase, to the same state that would 
be reached after a slow (adiabatic) process? This question has been answered recently 
in the affirmative for cooling expansions within the framework of the linear Schrodinger 
equation [5]. For preliminary work see 0, [7j. The method used to design the time- 
dependence of the trap frequency was based on Lewis-Riesenfeld invariants of motion 
[8] and simple inverse scattering techniques that had been applied for complex potential 
optimization [9J. Our objective here is to analyze if and how the same techniques used 
in that simple case can be adapted to non-linear interactions and systems described 
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by a Gross-Pitaevskii (GP) equation. As we shall see, the applicability of the method 
will depend critically on the effective dimension of the trap. We shall first discuss for 
simplicity with some detail one dimensional (ID) traps, and then 2D and 3D traps 
subjected to time-dependent frequencies. By ID traps we mean quasi-lD cigar-shaped 
traps with tight (fixed) transversal confinement where the axial frequency is varied 
in time; similarly 2D traps are quasi-2D disk-shaped traps with tight, fixed, axial 
confinement in which the transversal frequency is varied; and finally, the 3D traps 
refer to harmonic traps with spherical symmetry. We assume in all cases that a GP 
equation can be derived corresponding to each dimensionality, and use g generically for 
the coupling parameter of the non-linear term even though it is different for the three 
cases [TO] . 



2. One dimensional traps 



Our starting point is the effective ID Gross Pitaevski equation for the longitudinal (x) 
direction in an elongated cigar trap, 

dip 



ih 



dt 



■0. 



(i) 



g being the coupling parameter. The application of the invariant concept here is not as 
simple as for the Schrodinger equation [UJ, so we shall use instead an approach which 
leads in that case to the same results. The idea is to assume for the wavefunction the 
ansatz [121 



if)(x, t) = e 

Substituting this into Eq. 
$(p, t) = (f>(x,t), we get 
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where the dot means derivative with respect to time. Let us now impose that the 
coefficients in square brackets [...] of the last two terms vanish. This means that (we 
assume b real) 



a 1 l h imb 
p — - mo, a = — - 



2""~' " 2Kb 1 (4) 
and e ~( a + a *) x2 e -(P+P*) = -! Suppose now that the coefficient of 6 2 p 2 $ in ([3]) is 
made constant, equal to muj 2 /{2b A ) (for an alternative see the final discussion), where 
ujq = u>(0). Using PJ this is equivalent to imposing for b and u(t) an Ermakov equation, 



(5) 
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It is useful to express the resulting wave equation in terms of a new scaled time, 

and wavefunction *&(p, r) = (p(p,t), 

0¥ h 2 d 2 ^ muo 2 2t l|t|2t 

For g = this is the Schrodinger equation of a time-independent harmonic oscillator. 
The evolution of ip has thus been conveniently mapped to the simple solution of an 
auxiliary stationary system. Choosing 6(0) = 1, 6(0) = the "auxiliary" ([7]) and 
physical (CQ) oscillators coincide at t — 0, so any instantaneous eigenstate of t — 0, with 
vibrational quantum number n and energy E n = huoo(n + 1/2), evolves according to a 
propagating mode determined by Eqs. ( 12141) and the solution of the Ermakov equation 

m, 

if>(x, t) = 6- 1 / 2 e^^ 2 e-^ rW/ ^ n (x/6, 0). (8) 

In general this mode will not coincide with the instantaneous eigenstate of the physical 
Hamiltonian H(t) = — t^-§^ + ^muo(t) 2 x 2 , unless b(t) = [uoq/uo^)} 1 / 2 and b(t) = 0, up to 
the global phase factor e ~ lE ^ r ( t )/ n _ This motivated our proposal in [5j: it is an inverse 
method in which, given the initial uoq and final frequencies uuf = uo(tf), the intermediate 
trajectory u{t) is left undetermined at first and the boundary conditions 

6(0) = 1, 6(0) = 0, (9) 

Ktf) = hM) 1/2 , 6(f,) = (10) 

are imposed at initial and final times t = 0,t/ (they also imply a vanishing 6 at these 
two times to satisfy Eq. ([5])). b(t) is then interpolated with some functional form, 
e.g., a polynomial with enough coefficients to satisfy all conditions, and finally w 2 (t) 
is calculated from the Ermakov equation (J5j). This generates, in particular, very fast 
phase-space conserving cooling processes where u> 2 (t) takes during some time interval 
negative values, i.e., the trap becomes an expulsive potential. 

If g ^ the coefficient of the non-linear term in the auxiliary equation is generally 
time dependent. Thus, imposing b(tf) = eliminates the phase-factor e'^f^ 2 but 
nothing guarantees that ^/(r(tj)) is proportional to the instantaneous eigenstate of the 
GP equation at tf. A way out, in principle, is to make the coupling coefficient time- 
dependent with the aid of a Feshbach resonance as g(t) = go/b(t), with g constant. 
The resulting auxiliary equation has then time-independent coefficients, 

.8^ h 2 d 2 ^ muo 2 2t iti2t 

and can be solved in the form e~ l ^ T ^^ h ^/(x/b, 0), where \i is the chemical potential for 
the initial trap, so that 

il)(x, t) = 6- 1 /V^ x V^ rW/? ^(x/6, 0), (12) 
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and the same inverse method described for the Schrodinger equation can now be applied 
to design a fast frictionless process for the ground state condensate. One can easily 
check that, keeping b(t) = bf constant for t > tf, which results in u(i) = LOf and 
g = go^z/cjo) 1 / 2 for t > tf, the solution ip(x,t) of (pQ) given by (|T2j) becomes stationary, 
with a new scaled chemical potential fi/b(tf) 2 . 

Other special case is a "Thomas Fermi" (TF) limit, keeping g constant. Using a 
modified Ermakov equation and a different time scaling 

b + u(t) 2 b = ^, r(*) = jfy. (13) 

render an auxiliary equation with time-independent coefficients for the non-linear and 
harmonic potential terms. If g\fy\ 2 /(huJo) ^> 1 the kinetic term may be neglected, 

^ = m^ 

<9r 2 1 1 v ' 

This equation can also be solved by separation of variables, ^{x/b, r) = e~ l ^ T ^ h ^(x/b, 0), 
and ip(x,t) takes again the form of Eq. (Tl2l . with different values for /i, r, b, and the 
initial wavefunction. Note that this TF approximation is carried out in the auxiliary 
equation, and not at the level of the original GP equation, since that would imply 
a frozen density [321 [2]. From the modified Ermakov equation in ( ITBT) . the inversion 
method to find a frictionless trajectory u(t) requires in this 1D-TF scenario to change 
the boundary condition at tf in (flOl) to b(tf) = (cuo/ujf) 2 ^ 3 , with 6(0) = b(tf) = as 
before. 



3. Two and three dimensional traps 

The manipulations in ID suggest for 2D and 3D a wavefunction ansatz of the form [2] 
V>(r,t) = r d/2 e^5 0( r ,t), (15) 

where d is the dimension, r = (x 2 + y 2 ) 1 ^ 2 in 2D or r = (x 2 + y 2 + z 2 ) 1 ^ 2 in 3D. This 
form guarantees an auxiliary equation without first spatial derivatives. 

With p = r/b and a notation for the wavefunctions parallel to the ID case there 
results, by substituting fTl5|) into the 2D or 3D GP equations, 



%h— —b 2 = + — 

dr \dt 2m p 2 



p V* + -|-|*| 2 *, (16) 



b d- 

where the Laplacian should be adapted to the dimension, \1/ = ^f(p, r), and r has not 
been specified yet. (This equation includes the case d — 1 too by substituting r — > x 
and the Laplacian by a second derivative.) 

In 2D, the ordinary Ermakov equation (jSJ) and the r in Eq. ([6]) are the optimal 
choice since all coefficients in the auxiliary equation (assuming a constant g) become 
time independent, even outside the Thomas- Fermi regime, 

..dty h 2 . T mu 2 2t IT|2t 
or 2m 2 
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Figure 1. (color online). The squared frequency Lo 2 (t) for an expansion from 
uio = 250 x 2ir Hz to ujf = 2.5 x 2tt Hz in tf = 6 ms (a) Polynomial form b = Ylj=o a i^ ; 
(b) Exponential of a polynomial b — exp X)j=o c j'* J - both figures: ID, TF (solid, 
red line); 2D, or ID with g(t) = go/b(t), or 3D with g(t) = gob(t) (dotted, blue line); 
3D, TF (dot-dashed, magenta line). 



This is then the ideal situation for designing a frictionless process by shaping b and uj 
exactly as in the ID Schrodinger equation, i.e., using ([9]) and (jTUl) . 

Finally, the case d — 3 is considered. It is somewhat similar to ID, in the sense 
that the generic case leads to time- dependent coefficients in the auxiliary equation. 
Similarly to ID, by using Eqs. ([5]) and ([6]) the time-independence of the coefficients 
in the auxiliary equation would require now a time dependent coupling of the form 
9 if) = 9oKt)'i alternatively, in the Thomas-Fermi regime and with g constant, all 
coefficients become time independent with 

b + u{tfb=% r(t)= (18) 



b A ' w Jo b 3 

and in this case the boundary condition for b(tf) in fllUI) should be modified to 
b(tf) = (cuo/cuf) 2 / 5 , assuming again 6(0) = b{tf) = 0. 



4. Examples 




Figure 2. (color online), b corresponding to Fig. 1, b(tf) = (wo/w/) 2 ^", with v = 3 
(solid, red), v = 4 (dotted, blue), and v = 5 (dot-dashed, magenta), (a) Polynomial b, 
(b) b is the exponential of a polynomial. 
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Let us consider an expansion reducing the frequency 100 times from 250 x 2ir 
Hz to 2.50 x 2tt Hz in 6 ms. This time is too short for the condensate to follow 
any u(t) adiabatically at all times [5] but with our designed trajectories and thanks 
to the expulsive interval which accelerates the spreading, the final state would indeed 
be the same, up to a global phase, than the state obtained if a slow process could be 
implemented for such an expansion. For the regular Ermakov equation (jSJ), a polynomial 
b(t), and uq ^> uif, a simple estimate is that a repulsive time interval is necessary for 
tf < \/{2ujf). Figure 1 shows frequency trajectories for the different cases discussed 
before and two different ansatz of h. Higher powers of h in the right hand side of 
the Ermakov-like equations (corresponding to higher dimensions in the TF regime) 
imply a smaller increment for b during the expansion, see Fig. 2, which makes the 
change of u smoother as well. It is remarkable that for a fixed g (for 2D, or the 
TF regimes in ID and 3D), the non-linearity does not play any role in the design of 
optimal (frictionless) frequency trajectories. They only depend on the initial and final 
frequencies, the available time tf and the functional form chosen for b(t). 



5. Discussion 



We will provide here some complementary details and relate the results to other works. 
An important remark on the TF approximation used for ID and 3D geometries is that 
the non-linear coupling cannot be arbitrarily strong. The condition g\ty\ 2 /(huJo) 1 
should be compatible with the derivation of the ID GP equation [10J in a weak 
interaction limit, i.e., a s |?/>| 2 << 1, where a s is the s-wave scattering length. 

For completeness we should mention an alternative to the steps given after Eq. (j3J). 
We may as well impose that the coefficient multiplying p 2 b A ^> must vanish instead of 
becoming a non-zero constant p31 [HI [15]. This amounts to imposing b + uj{t) 2 b = 
instead of the Ermakov equation (jSJ). Proceeding as in Sec. [2] with r given by Eq. ([6]), 
the resulting auxiliary equation becomes 

fl* h 2 d 2 ^ JIT|2t 

* n 87 = -2^W + 9m% (19) 
which is not an equation for the harmonic oscillator but for a condensate without 
confining external fields and with a, generically, time dependent non-linear coupling 
factor. Adapting the time dependence of g as g(t) = go/b(t), this method provides, 
from known analytical solutions of Eq. ( fl9|) with a constant factor g(t)b(t) = go, explicit 
solutions that have been used in the context of soliton dynamics [T3l [T^ [15]. While the 
solutions ip(x,t) for the same u(t) and initial conditions should of course be equivalent 
to the ones obtained with the ordinary Ermakov equation, we find the later better suited 
for the application of our inverse technique. 

In summary, it is possible to take a Bose-Einstein condensate in a very short time 
from an initial harmonic trap to a final one without excitations, by choosing the time 
dependence of the frequency according to the Ermakov equation or its modifications 
after matching the time dependence of a scaling factor to suitable boundary conditions. 
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In ID and 3D traps this requires either a simultaneous change of the time-dependence 
of the coupling, or a Thomas-Fermi type of regime. 2D traps are privileged in this 
respect and do not require either of these conditions. Their peculiar symmetry properties 
were already noticed by Pitaevskii and Rosch [16]. Indeed, the 2D geometry allows for 
an extension of the present results beyond the GP equation framework by expanding 
perturbatively the field operator around the condensate wavefunction, and treating the 
perturbation with an ansatz parallel to (fT5j) and the same phase [2]. 
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